home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / VideoToolboxSources / Normal.c < prev    next >
C/C++ Source or Header  |  1995-07-20  |  8KB  |  238 lines
  1. /*
  2. Normal.c
  3. statistical functions related to the normal distribution.
  4. Also see: Binomial.c,ChiSquare.c, Exponential.c, Uniform.c
  5.  
  6.     f=NormalPdf(x);
  7.     p=Normal(x);
  8.     x=InverseNormal(p);
  9.     x=NormalSample();
  10.     f=Normal2DPdf(double r);
  11.     p=Normal2D(r);
  12.     r=InverseNormal2D(p);
  13.     r=Normal2DSample();
  14. Normal2D is a Gaussian pdf over two dimensions, integrated over all orientations,
  15. 0 to 2π. r is the distance from the origin, [0,Inf]. 
  16.  
  17.     BoundedNormalIntegers(distribution,n,mean,sd,min,max);
  18. Fills the "distribution" array with n ordered integers so that random samples
  19. from the array,
  20.     i=distribution[nrand(n)];
  21. will have, as nearly as possible, the specified distribution, i.e. they will
  22. be samples (rounded to integer) drawn from the interval [min-.5,max+.5],
  23. where min and max are integers, of a normal distribution with the specified
  24. mean and variance. Specifying an infinite sd will result in a uniform
  25. distribution over the interval [min,max]. (This is detected as a special
  26. case and performed quickly.) Once the distribution array has been filled,
  27. random samples from it are a fast way of generating bounded Gaussian noise
  28. to be added to each pixel of an image. Note that "mean" and "sd" apply only
  29. to the unclipped underlying normal distribution. Call this
  30.     mean=MeanW(distribution,n,&sd);
  31. to compute the mean and sd of your integer distribution. The runtime of 
  32. BoundedNormalIntegers will be proportional to max-min+1, and nearly independent
  33. of n. It takes about 0.3 s on a Mac IIci for max-min+1=256.
  34.  
  35. Copyright (c) 1989,1990,1991,1992,1995 Denis G. Pelli
  36. HISTORY:
  37. 1989    dgp wrote it.
  38. 4/8/90    dgp    changed the names of the routines. 
  39.             Made sure that domain error produces NAN.
  40. 6/90    dgp    added NormalSample()
  41. 7/30/91    dgp    now use NAN defined in VideoToolbox.h
  42. 12/28/91 dgp sped up NormalPdf() by calculating the scale factor only once
  43. 12/29/91 dgp extracted code to create new routine UniformSample.c
  44. 1/11/92    dgp    rewrote Normal()'s polynomial evaluation to halve the number of multiplies
  45.             Renamed NormalPDF() to NormalPdf().
  46. 1/19/92    dgp    defined the constants LOG2 and LOGPI in VideoToolbox.h
  47.             Added more domain tests, returning NAN if outside. 
  48.             Added more checks to main().
  49.             Wrote Normal2DPdf(),Normal2D(),InverseNormal2D(),and Normal2DSample().
  50. 2/1/92    dgp    Redefined Normal2DPdf(r) to now, more sensibly, treat r as the random
  51.             variable, rather than the implicit x and y, where r=sqrt(x^2+y^2). 
  52.             Previously, Normal2D(R)=Integrate[2 Pi r Normal2DPdf(r),{r,0,R}]
  53.             now Normal2D(R)=Integrate[Normal2DPdf(r),{r,0,R}]. There is no change
  54.             in Normal2D(). I suspect that Normal2D is in fact the Raleigh distribution,
  55.             and I will rename it if this is in fact the case.
  56. 11/13/92 dgp InverseNormal(0) now returns -INF, and InverseNormal(1) returns INF.
  57. 12/15/93 dgp declared some arguments "register".
  58. 3/26/94    dgp    added BoundedNormalIntegers().
  59.             Asked compiler to use 68881 instructions for transcendentals.
  60. */
  61. #include "VideoToolbox.h"
  62. void BoundedNormalIntegers(short *distribution,long n,double mean,double sd
  63.     ,short min,short max);
  64.  
  65. double NormalPdf(register double x)
  66. /* Gaussian pdf. Zero mean and unit variance. */
  67. {
  68.     if(IsNan(x))return x;
  69.     return exp(-0.5*(x*x+(LOG2+LOGPI)));
  70. }
  71.  
  72. double Normal(register double x)
  73. /*
  74. Cumulative normal distribution. From Abramowitz and Stegun Eq. (26.2.17).
  75. Error |e|<7.5 10^-8
  76. */
  77. {
  78.     register double P,t;
  79.     
  80.     if(x<0.0) return 1.0-Normal(-x);
  81.     t=1.0/(1.0+0.2316419*x);
  82.     P=(0.319381530+(-0.356563782+(1.781477937+(-1.821255978+1.330274429*t)*t)*t)*t)*t;
  83.     return 1.0-NormalPdf(x)*P;
  84. }
  85.  
  86. double InverseNormal(register double p)
  87. /*
  88. Inverse of Normal(), based on Abramowitz and Stegun Eq. 26.2.23.
  89. Error |e|<4.5 10^-4.
  90. */
  91. {
  92.     register double t,x;
  93.     
  94.     if(IsNan(p))return p;
  95.     if(p<0.0)return NAN;
  96.     if(p==0.0)return -INF;
  97.     if(p>0.5) return -InverseNormal(1.0-p);
  98.     t=sqrt(-2.0*log(p));
  99.     x=t-(2.515517+(0.802853+0.010328*t)*t)/(1.0+(1.432788+(0.189269+0.001308*t)*t)*t);
  100.     return -x;
  101. }
  102.  
  103. double NormalSample(void)
  104. {
  105.     return InverseNormal(UniformSample());
  106. }
  107.  
  108. double Normal2DPdf(double r)
  109. /* Gaussian pdf over two dimensions, integrated over all orientations, 0 to 2π. */
  110. /* The argument is taken to be the distance from the origin, [0,Inf]. */
  111. /* The rms is 1 */
  112. {
  113.     if(IsNan(r))return r;
  114.     if(r<=0.0)return 0.0;
  115.     return 2*r*exp(-r*r);
  116. }
  117.  
  118. double Normal2D(double r)
  119. /* Integral of Normal2DPdf() from zero to r. */
  120. {
  121.     if(IsNan(r))return r;
  122.     if(r<=0.0)return 0.0;
  123.     return 1.0-exp(-r*r);
  124. }
  125.  
  126. double InverseNormal2D(double p)
  127. {
  128.     if(IsNan(p))return p;
  129.     if(p<0.0 || p>1.0)return NAN;
  130.     return sqrt(-log(1.0-p));
  131. }
  132.  
  133. double Normal2DSample(void)
  134. /* rms is 1 */
  135. {
  136.     return InverseNormal2D(UniformSample());
  137. }
  138.  
  139. /*
  140.     Integrate[exp(-.5*u^2),{u,a,a+1/sd}]
  141.     =exp(-.5*a^2)*Integrate[exp(-.5*((a+e)^2-a^2)),{e,0,1/sd}]
  142.     =exp(-.5*a^2)*Integrate[exp(-.5*e*e)*exp(-a*e),{e,0,1/sd}]
  143.     =exp(-.5*a^2)*Integrate[(1-.5*e*e)*exp(-a*e),{e,0,1/sd}]
  144.     =exp(-.5*a^2)*((exp(-a/sd) - 1)/(-a)-.5*Integrate[e*e*exp(-a*e),{e,0,1/sd}])
  145.     =exp(-.5*a^2)*(1-exp(-a/sd))/a
  146. */
  147. void BoundedNormalIntegers(register short *distribution,long n,double mean,double sd
  148.     ,short min,short max)
  149. {
  150.     register short i;
  151.     register long j,count,values,round;
  152.     double p,p0,p1,x;
  153.     short shortcut;
  154.     
  155.     j=0;
  156.     if(IsInf(sd)){
  157.         /* Uniform distribution over the interval [min,max] */
  158.         values=max-min+1;
  159.         assert(n<LONG_MAX/values);
  160.         round=values/2;
  161.         for(i=min;i<max;i++){
  162.             count=((i-min+1)*n+round)/values;
  163.             while(j<count)distribution[j++]=i;
  164.         }
  165.     }else{    
  166.         shortcut=sd*sd>n;    /* guarantees count will err by at most ±0.5 */
  167.         p=p0=Normal((min-.5-mean)/sd);
  168.         p1=Normal((max+.5-mean)/sd)-p0;
  169.         for(i=min;i<max;i++){
  170.             x=(i+.5-mean)/sd;
  171.             if(shortcut)p+=NormalPdf(x)*(exp(x/sd)-1)/x;
  172.             else p=Normal(x);
  173.             count=0.5+n*(p-p0)/p1;
  174.             while(j<count)distribution[j++]=i;
  175.         }
  176.     }
  177.     while(j<n)distribution[j++]=max;
  178. }
  179.  
  180.  
  181. #if 0 /* A test program. */
  182.     void main()
  183.     {
  184.         double x,y,sum,dx,a,b,mean,sd;
  185.         static double z[1000];
  186.         int i;
  187.         
  188.         Require(0);
  189.         srand(clock());
  190.         printf("%4s%15s%15s%20s%15s\n","x","NormalPdf(x)","Normal(x)","InverseNormal","Error");
  191.         for(x=-4.0;x<=4.0;x+=2.0){
  192.             printf("%4.1f%15.8f%15.8f%20.4f%15.4f\n",
  193.             x,NormalPdf(x),Normal(x),InverseNormal(Normal(x)),InverseNormal(Normal(x))-x);
  194.         }
  195.         sum=0.0;
  196.         dx=0.001;
  197.         for(x=-1.;x<0.;x+=dx)sum+=NormalPdf(x);
  198.         sum*=dx;
  199.         sum-=Normal(0.0)-Normal(-1.0);
  200.         printf("Partial integral of NormalPdf error %.5f\n",sum);
  201.         for(i=0;i<1000;i++)z[i]=NormalSample();
  202.         mean=Mean(z,1000,&sd);
  203.         printf("1000 samples mean %.2f sd %.2f\n",mean,sd);
  204.         printf("\n");
  205.     
  206.         printf("%4s%15s%15s%20s%15s\n","x","Normal2DPdf(x)","Normal2D(x)","InverseNormal2D","Error");
  207.         for(x=-1.;x<=5.0;x+=1.0){
  208.             printf("%4.1f%15.8f%15.8f%20.4f%15.4f\n",
  209.             x,Normal2DPdf(x),Normal2D(x),InverseNormal2D(Normal2D(x)),InverseNormal2D(Normal2D(x))-x);
  210.         }
  211.         sum=0.0;
  212.         dx=0.0001;
  213.         for(x=0;x<1.;x+=dx)sum+=Normal2DPdf(x);
  214.         sum*=dx;
  215.         sum-=Normal2D(1.0);
  216.         printf("Partial integral of Normal2DPdf error %.5f\n",sum);
  217.         for(i=0;i<1000;i++)z[i]=Normal2DSample();
  218.         mean=Mean(z,1000,&sd);
  219.         printf("1000 samples rms %.2f\n",sqrt(mean*mean+sd*sd));
  220.         printf("\n");
  221.         for(i=0;i<1000;i++){
  222.             x=NormalSample();
  223.             y=NormalSample();
  224.             z[i]=sqrt((x*x+y*y)/2.);
  225.         }
  226.         mean=Mean(z,1000,&sd);
  227.         printf("1000 (x,y) normal samples with sd 2^-0.5 have rms hypotenuse of %.2f\n",sqrt(mean*mean+sd*sd));
  228.         printf("\n");
  229.     
  230.         a=4.0*atan(1.0);
  231.         if(a!=PI)printf("4*atan(1)-PI %.19f\n",a-PI);
  232.         a=log(a);
  233.         if(a!=LOGPI)printf("Error: log(PI) %.19f, error in LOGPI %.19f\n",a,LOGPI-a);
  234.         a=log(2.0);
  235.         if(a!=LOG2)printf("Error: log(2) %.19f, error in LOG2 %.19f\n",a,LOG2-a);
  236.     }
  237. #endif
  238.